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Abstract. Many schemes have been proposed to perform a model-independent constraint 
on cosmological dynamics, such as nonparametric dark energy equation of state (EoS) uj(z) 
or the deceleration parameter q(z). These methods usually contain derivative processes with 
respect to observational data with noise. However, it still remains remarkably uncertain 
when one estimates the numerical differentiation, especially the corresponding truncation 
errors. In this work, we introduce a global numerical differentiation method, first formulated 
by Reinsch(1967), which is smoothed by cubic spline functions. The optimal solution is 
obtained by minimizing the functional $(/). To investigate the potential of the algorithm 
further, we apply it to the estimation of the transition redshift zt with simulated expansion 
rate E{z) based on observational Hubble parameter data(OHD). An effective method to 
determine the free parameter S appearing in Reinsch Splines is provided. 
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1 Introduction 

The central task of modern cosmology is to uncover the dynamic evolution and the geometric 
structure of the universe. According to plenty of cosmological observations, such as distant la 
supernovae (SNe la), cosmic microwave background (CMB) and so forth, the recent universe, 
dominated by the so-called dark energy, is undergoing accelerated expansion. The Lambda 
CDM model (ACDM) based on general relativity can provide consistent explanation with 
the observational data. Whereas, the nature of the dark energy still remains mysterious. In 
order to study the characteristics of the accelerated universe, many different methods have 
been proposed to free the limitation of the concrete cosmological models, e.g. the direct 
reconstruction of the equation of state (EoS) of dark energy uj(z) [1, 2] or the expansion rate 
and deceleration parameter E(z), q(z) [3]. 

How to treat the numerical differentiation correctly is one common problem appearing 
in the above works. The accuracy of numerical differentiation of noisy observational data is 
difficult to control. It encompasses many subtleties and pitfalls that may cause large error in 
the actual computation. From the mathematical point of view [4], the derivative of a given 
function is obtained by infinitesimal calculus, it's, however, impractical for real experimental 
data due to the discrete property. Additionally, for most cases, we must take the measured 
error or noise into consideration. An effective estimation of the derivative of noisy sample 
y(xi) with respect to the variable xi has to cope with these two restrictions. 

To compute the numerical differentiation of a given sample without knowing the under- 
lying function, typically one needs to obtain the approximation with some basic functions. 
Then one can hope that the derivative of the approximation is good enough to represent the 
actual situation. Finite differences is the general, but very crude method which works by 
fitting or interpolating on some sub-interval of the sample domain. If the boundary points 
of the sub-interval are very close, the subtraction of the two numbers will be ill-conditioned. 
It may render the final result meaningless. In addition, the optimal function obtained by 
such method is not smooth. As a result, it's hard to represent the real world. Therefore, 
finite differences method may be a bad choice in many cases. Reinsch [5] developed another 
optimal algorithm to perform numerical differentiation with spline functions, here referred as 
Reinsch Splines. The method can overcome some difficulties arising from the Finite differ- 
ences. The main task of the work is to analysis the Reinsch Splines and its practical potential 
on observational data. More detailed descriptions will be presented in the following sections. 
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This paper is organized as follows: The fundamental algorithm of Reinsch Splines and 
the corresponding errors are furnished in Sec. 2. In order to investigate the potential of 
the algorithm, we apply the algorithm in Sec. 3 to estimate the transition redshift Zt with 
simulated cosmological expansion rate E{z) based on OHD. Finally, the brief summary and 
discussion will be given in Sec. 4. 

2 Reinsch Splines Method 

Reinsch Splines is first proposed by Reinsch [5] to replace strict interpolation by some kind 
of smoothing. The appropriate trial functions to estimate the experimental data is spline 
functions. 

Given a sample (xi, (i = 1, ...,n) which satisfies that 

Xl < x 2 < ... < x n , 

then the problem can be stated as a special instance of Tikhonov regularization method, 
which is just to look for the minimum of the functional 

*(/) = « |e ( Vt ~l (Xl) ) 2 - + limn 2 , (2.i) 

where Oi is the noise, f(x) is the optimal square integrable function over the domain. | \f"(x)\ | 
denotes the L 2 -norm 



\\f"(x)\\ = (£ n f"(x) 2 dx 



1/2 



S > is a given constant, allowing for an implicit rescaling of the quantities <7j, which controls 
the extent of smoothing. If <t, is the estimate of the standard deviation of yi, the value will 
lie within 

n - V2n < S <n + V2n. (2.2) 

It should be noted that S is very sensitive to the experimental noise crj, sometimes the 
limitation must be relaxed, especially when the noise is crude, although the confidence interval 
of S is accurate enough for small deviations, a is the Lagrangian parameter satisfying 
d<S>(f)/da = and a / 0. 

Hanke & Scherzer [6] provides a rigorously proof that the minimizer of equation (2.1) 
is a natural cubic spline. Following the notation of Reinsch [5], we express f(x) as 

f(x) = ai + bi(x - Xi) + Ci(x- Xi) 2 + di(x- Xi) 3 , (2.3) 

where 

Xi < x < Xi + \ , i = 1, n — 1. 

The spline functions satisfy specified smoothing and boundary conditions. But there is no 
need to agree exactly with the experimental data yi on the node Xj. Once we obtain the 
coefficients aj, bi, Ci and di, the spline functions will be determined uniquely. A constructive 
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algorithm for calculating the splines has been given by Reinsch [5]. In the remaining part of 
this section, we will focus on the error analysis of the Reinsch Splines. 

Three sources of errors will impact the final evaluation of the numerical differentiation: 
experimental errors resulting from individual data point, truncation errors between the op- 
timal f'(x) and the true derivatives, and the rounding errors due to the narrow length of 
sub-interval. It is precisely because of the existence of the rounding errors attributed to the 
loss of significance of digits, finite differences fails to estimate the differentiation effectively. 
Reinsch Splines, however, can overcome the embarrassment to some extent. The truncation 
errors greatly dependent on the choice of the algorithm and the fitting accuracy. Hanke &; 
Scherzer [6] provides a formula of the truncation errors, and the rigorous proof can also be 
found. Note that some slight adjustments should be done to match the Reinsch Splines, 
though the correction is straightforward. As for the experimental errors, the propagation of 
them is dealt with the conventional expression 

* <2 - 4) 

where m denotes the coefficients of the spline functions, and iV is the number of the points 
of sub-interval. 



3 Estimation of the Transition Redshift 



Observations show that the universe has undergone a dynamic phase transition from de- 
celeration to acceleration expansion which leads to the change of the sign of deceleration 
rate q(z). The transition redshift z% has profound impact on the evolution of the universe. 
Different groups have measured Zt with SN la , BAO and CMB observations [7, 8]. In or- 
der to constrain zt model-independently, various parametric expressions of q(z) have been 
presented. Such parametric expressions e.g. oj{z) and q(z), are extremely convenient and 
effective for cosmological research. But the situation will remarkably change if one introduces 
the nonparametric forms of them. Just as we have mentioned, derivative calculations often 
appear in these expressions. Therefore the first difficulty must to be solved is the accurate 
estimation of the numerical differentiation with noise observational sample. 

Recently, observational Hubble parameter data (OHD) H (z) has obtained many atten- 
tions on the constraint of cosmological parameters [9-11]. Ma & Zhang [12] and Zhang et 
al. [13] have summarized the power and potential of OHD from statistical point of view. 
At present, there are total about 28 independent measurements of H(z), which are listed in 
table 1. 

The relation between transition redshift and H(z) can be written as 



z t 



1 dH(z) 
H(z) dz 



1. (3.1) 



z=z t 



Obviously, the technique of numerical differentiation needs to be used to determine Zt if 
H(z) sample is available. Unfortunately, unlike the SN la dataset, the recent OHD are too 
penurious to provide valuable estimation of the transition redshift using such method. There 
are some other reasons that make it invalid. The apparent uncertainty of the existing H{z) 
will greatly increase the possibility of failure. Moreover, sparse sample will lead the optimal 
functions to deviate the actuality seriously, as well as the corresponding differentiation. Lima 
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Figure 1. Top: offsets of E der (z) with r/ = 0.1320. The solid line is \e(z)\ = r)E fid (z). Most of the 
derived--E(z) follow the relation |e| < r]E der (z). Bottom: mean value boundary of E(z). Two dashed 
curves just cover the very most E der (z) points of interest. Obviously, most of the derived- E(z) locate 
within Efid + £- < E der < E fld + e + . 

et al. [14] suggests three different means to obtain H(z), so it is expected that such embar- 
rassments could be overcome by ongoing and future observations to expand the volume of 
the sample. 

3.1 Simulated E(z) 

To extend the applications of OHD further and test the Reinsch Splines, a simulated sample 
of H{z) will be helpful. Further, from equation (3.1) the dimensionless expansion rate E(z), 
defined as H(z)/Hq, will be a better choice for the present analysis. Therefore, we will 
concentrate on the technique of generating a sample of E(z), named as simulated- in 
this section. Our simulation is based on the spatially flat ACDM model with S} m = 0.28 and 
SIa = 0.72. A Gaussian prior Hq = 74.2 ± 3.6 km-s~ 1 -Mpc~ 1 suggested by Riess et al. [15] is 
adopted. The expansion rate E(z) in the fiducial model can be written as 

E fid (z) = y/Q m (l + zf + n A . (3.2) 

Based on the OHD listed in table 1 and equation (2.4), we can obtain the 'observational' 
expansion rate, denoted as derived- E(z), and the corresponding error. 
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Figure 2. The model of <J S i m (z). The circles represent the derived- E^z) sample where the solid ones 
are chosen to fit the up- and down-edge. The solid line is the midline of the two edges. 

The simulated-i?^) sample should take the fiducial model value Et^ as its expectation 
and follow the same systematic information and characteristic of the derived- E(z). However, 
this goal is so difficult to reach, because we have quite a few real data points and almost 
little knowledge about the derived- E(z), e.g. the distribution of the data points along the 
redshift-axis. Therefore, our E(z) simulation includes the following two procedures: 

The offset estimation e(z): To generate the simulated- E(z) sample at any given redshift 
point, we introduce the variable e(z) satisfying 

E sim (z) = E fld (z) + e(z), (3.3) 

which represents the offset between the fiducial value and simulation at given z , In order to 
make (E S i m ) = Ef i( i, e(z) should be a random variable with respect to z. 

Defining the offset i(z) = E c i er (z) — Efid(z) whose absolute value is shown in figure 
1. Assuming that there is no bias on the sign of it, then we can make the offset satisfy 
\i(z)\ < E(z)rj for most derived- E(z) points, where r\ is a constant. The best estimate value 
of T] is 0.1320. From figure 1 we see that the offset curve covers most of the points of interest, 
ignoring only three points with extraordinary deviation. As a result, we can denote 

e±(z) = ±E fid (z)r t ( V = 0.1320) 
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Table 1. The currently known observational Hubble parameter data. The latest four data points 
from Zhang et al. [16] are added. However, we do not include them to perform the simulation because 
of the remarkable uncertainty and the existing dense data points with relatively small deviations at 

low redshift. 

redshift z H(z) &h(z) 

(km • s^ 1 • Mpc" 1 ) (km • s _1 • Mpc~ l ) 



0.090 69 12[17] 

0.170 83 8[17] 

0.270 77 14[17] 

0.400 95 17[17] 

0.900 117 23[17] 

1.300 168 17[17] 

1.430 177 18[17] 

1.530 140 14[17] 

1.750 202 40[17] 

0.480 97 62[18] 

0.880 90 40[18] 

0.179 75 4[19] 

0.199 75 5[19] 

0.352 83 14[19] 

0.593 104 13[19] 

0.680 92 8[19] 

0.781 105 12[19] 

0.875 125 17[19] 

1.037 154 20[19] 

0.24 79.69 3.32[20] 

0.43 86.45 3.27[20] 

0.44 82.6 7.8[21] 

0.60 87.9 6.1[21] 

0.73 97.3 7.0[21] 

0.07 69.0 19.6[16] 

0.12 68.6 26.2[16] 

0.20 72.9 29.6[16] 

0.28 88^ 36.6[16] 



as the boundaries of the offset values of simulated expansion rate (see figure 1). As a random 
variable, e(z) follows the Gaussian distribution N(0, r qEfi ( [(z)/2) so that the probability of 
e(z) falling within the offset domain is 95.4%. 

The uncertainty estimation <j(z): There is an apparent trend of the errors of derived- 
E(z): the uncertainty 0e(z) becomes larger as the increase of redshift z except for two outliers 
at z = 0.48 and 0.88, as figure 2 shows. Ignoring the two points, we find the rest are right 
within the region between the lines 

<t+ = 0.2324z + 0.1365, cr_ = 0.1091z + 0.0393. 
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And the midline between the two boundaries is 



0", 



2^+ + CT -) 



which represents the expectation value of the uncertainty cr s i m of the simulated expansion 
rate. Meanwhile, we make the a s i m follow a Gaussian distribution N (a m (z) , p(z)) at any 
given redshift, where 



ensuring that the simulated uncertainty a s im falls within the region between lines cr+(z) and 
a~(z) with the probability of 95.4%. 

According to the above two procedures, the method of generating a simulated E(z) 
sample is completed. First of all, we can calculate the fiducial values Et^z) from equation 
(3.2). A random variable e{z) can be drawn from a Gaussian distribution, so the E s i m (z) is 
generated via equation (3.3) at a given z. Finally, the corresponding uncertainty cr s i m (z) is 
also estimated through another Gaussian distribution. Some similar but different techniques 
have also been proposed by Ma & Zhang [12] and Wang & Zhang [22] to simulate OHD. Figure 
3 shows the final simulation which contains total 256 data points ranging from 0.001 < z < 
2.0, as well as 24 derived--E(z). 

3.2 Numerical Results 

In this section, the Reinsch Splines is applied to estimate the transition redshift with the 
simulated expansion rate E(z). 

In order to obtain the optimal function f(x) defined as equation (2.3), we need to 
minimize the functional <&(/) (see equation (2.1)). For such noisy data, constant S plays an 
important role. A proper choice of the constant will substantially boost the fitting accuracy 
of numerical differentiation. If the underlying function to be fitted is known, the classical 
least-squares fitting can be used. But how to determine it without any model information, 
just as the model-independent tests in cosmology, and how to judge whether the chosen S 
can approach the real situation remain to be solved. The general practice is to set S — n 
which satisfies equation (2.2). It works fairly well when we test the method on different 
analytic functions with small errors. An example can also be seen Reinsch [5]. But it fails 
for our simulated sample because of the significant errors of E[z)\ 

One possible estimation of S here we introduce is to minimize 



where Vj = (xj, f'(xi)). With this condition, we can make the numerical differentiation as 
smooth as possible. However, the criterion can merely ensure local smoothing. In other 
words, the global oscillations or undulations due to the accuracy of Reinsch Splines can not 
be determined by equation (3.4). In most cases, we hope that the number of times such 
oscillations or undulations happen should be as little as possible. One convenient way is to 
calculate the number of the extremum values iV e when different S is given. We can find N e 
using the general definition of extremum on mathematics. Equivalently, it can be determined 



n-l 




(3.4) 



i=i 
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Figure 3. Simulated E(z) sample based on fiducial ACDM model. The 24 OHD derived- E(z) points 
are also displayed with error bar. The solid line represents the E(z) in ACDM with VL m = 0.28 and 
n A = 0.72. 

through the 'nodes', corresponding to the locations of extremum or inflection points, in the 
region of errors, just as figure 4 shows. Utilizing the two criterions, the constant S can be 
evaluated uniquely. In this work, the best value is S = 76.97. 

Once S is fixed, the optimal f(x) can also be obtained with the algorithm presented in 
section 2. Figure 4 illustrates the final numerical results. The best estimate of the transition 
redshift is z% = 0.69^o'?4 which is moderately consistent with the observational and theoretical 
results. Note that we just take the statistical errors and truncation errors into consideration 
due to the rounding errors are negligible in our sample. Furthermore, from figure 4, we see 
that there are large deviations of the boundary ascribed to the selection of the boundary 
conditions in Reinsch Splines. However, a different version of the boundary conditions based 
on physical facts or some prior knowledge may be helpful to overcome the problem. 

4 Conclusions 

In the present work, we introduce a general numerical differentiation referred as Reinsch 
Splines, and preliminarily apply it to the analysis of the transition redshift zt with simulated 
expansion rate E(z). By comparing the difference between the model value of Zt and the 
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Figure 4. The determination of the transition redshift using simulated E(z). The error bar contains 
both statistical errors and truncation errors. The transition redshift z t ~ 0.69 estimated using Reinsch 
Splines. The 'nodes' showing in the region of errors correspond to the locations of extremum or 
inflection points of the differentiation. 

estimation, some experimental characteristics of the Reinsch Splines are unfolded. The deter- 
mination of constant S is essential to approach the optimal f(x), and we provide an effective 
recipe to evaluate it when the experimental error a is significant. Meanwhile, Reinsch Splines 
can also be applied to reconstruct the nonparametric dark energy equation of state uj(z) with 
recent supernova data. Some different techniques have been provided by Clarkson & Zunckel 
[1] and Holsclaw et al. [2]. By contrast, Reinsch Splines will generalize the procedures of 
numerical differentiation and reduce the artificial uncertainty. 

The accuracy of the Reinsch Splines still remains to improve to extend its applications. 
One feasible method is to relax the boundary conditions when constructing the optimal spline 
functions f(x). In addition, the errors resulting from the numerical differentiation should 
be dealt with carefully, and need further discussion, especially the truncation errors and 
rounding errors. 
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